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Abstract 

The stochastic differential equations for a model of dissipative particle 
dynamics, with both total energy and total momentum conservation at every 
time-step, are presented. The algorithm satisfies detailed balance as well as the 
fluctuation- dissipation theorems that ensure that the proper thermodynamic 
equilibrium can be reached. Macroscopic equilibrium probability distributions 
as well as equations of state for the model are also derived, and an appropriate 
definition of the free energy of the system is consistently proposed. Several 
simulations results of equilibrium as well as transport properties, including 
heat transport and thermal convection in a box, are shown as proof of the 
internal consistency of the model. 
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1 Introduction 



The computer simulation strategy for dynamics of complex systems known as Dis- 
sipative Particle Dynamics or, simply, DPD, has been the subject of several studies 
in recent years. This methodology was introduced by Hoogerbrugge et a/.|l[ to 
model the dynamic behaviour of fluids by using a particulate method in which an 
ensemble of mesoscopic particles interact with each other via conservative as well as 
dissipative forces. Contrary to other particulate Lagrangian met hods 0], oriented at 
modelling the hydrodynamics of macroscopic systems, DPD also incorporates Brow- 
nian forces in the particle-particle interactions. In this way, thermal fluctuations can 
be described and a certain thermodynamic equilibrium exists. DPD is, therefore, 
especially suited to the modelling of mesoscopic systems, where fluctuations play an 
important role. One can mention, for instance, the dynamics of polymer solutions, 
nucleation and phase separation problems, dynamics of confined fluids, colloidal sus- 
pensions, and kinetics of chemical reactions in confined geometries, etc. The DPD 
method has already been succesfully applied to some cases of practical interest such 
as to the rheology of colloidal suspensions 0] and to the microphase separation of 
diblock copolymer solutions 0, among others. 

In DPD, dissipative and random interactions are pairwise and chosen in such a 
way that the center of mass motion of each interacting pair is insensitive to these fric- 
tional and random forces. Hence, on the one hand, if closed and thermally isolated, 
the system relaxes fast to its thermal equilibrium and, on the other, its overall mo- 
mentum is a conserved variable in view of the particle's momentum exchange rules. 
This second feature allows the system to exhibit a hydrodynamic behaviour from 
a macroscopic point of view, that is, for length and time scales larger than those 
characterizing the particle-particle interactions. With respect to other mesoscopic 
models for hydrodynamic behaviour such as Lattice Gas[|3| or Cellular Automata[|J, 
the DPD model is isotropic and Galilean-invariant due to the fact that it is not de- 
fined on a lattice and is straightforwardly based on Newton's equations of motion, 
as in MD. In addition, it has no extra conservation laws and is also computationally 

2 



efficient. As originally formulated, however, the DPD model can only deal with 
isothermal conditions since dissipative and Brownian forces cause no energy conser- 
vation in the particle-particle interaction. The equation for the energy transport is 
then of the relaxation type|f§ and, therefore, heat flow, related to the energy con- 
servation at the microscopic scale, cannot be described within the framework of the 
original model (which will be referred to as isothermal DPD from now on). In many 
problems of interest, either fundamental or applied, the study of the transport prop- 
erties of heat at the mesoscopic level are very important. Thus, the incorporation 
of the thermal effects into a DPD algorithm is necessary for future applications of 
the method to problems of relevance. 



In a previous paper |T]| we introduced an algorithm in which the conservation 



of the total energy in the particle-particle interaction was consistently taken into 
account with the inclusion of the particle's internal energy in the model (this model 
will be referred to as DPDE from now on). In this paper we will analyze several 
aspects of this novel DPDE algorithm. In the first place, we discuss the extension 
of our original algorithm, given in ref.]18[|, to incorporate, on the one hand, temper- 



ature dependencies in the dynamic parameters of the model so that the treatment 
of arbitrary temperature-dependencies in the transport coefficients can be treated, 
something lacking in the older DPD models. On the other hand, we develop a direct 
derivation of the algorithm from the Langevin equations of motion for the relavant 
variables, by introducing a interpretation rule for these Langevin equations that 
ensures that the detailed balance condition is satisfied in an Euler-like algorithm. 
The detailed balance condition is necessary for the system to evolve towards the 
proper thermodynamic equilibrium. Furthermore, the algorithm obtained in this 
way conserves the energy at every time-step instead of in the mean as was the case 
in the older DPDE algorithm |18||. We have noticed that different algorithms lead 



to the same Fokker-Planck equation, i.e., different stochastic processes can lead, in 
fact, to the same dynamics for the probability distribution. In the second place, 
we have analyzed the macroscopic equilibrium properties of the DPDE system. We 
have shown that a thermodynamic analysis of the DPDE system is possible in view 
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of the fact that the equilibrium probability distribution is known. Properties such 
as a free energy or the equations of state can be derived in terms of the parameters 
defining the model. The simulation results of the equilibrium properties show an 
excellent agreement with the theoretical predictions. Third and last, we have done 
simulations in systems under a temperature gradient imposed by the temperature of 
two walls. We have found temperature profiles, as expected since the model allows 
the simulation of heat transport, and the existence of inhomogeneous temperature 
distributions, which is one of the main virtues of the DPDE algorithm. Furthermore, 
if a gravity field is considered, the system undergoes convective motion in addition 
to the temperature gradient, the Rayleigh number qualitatively estimated as being 
of the order of 10 3 . 

The paper is organized as follows. In section II, we introduce the main hypoth- 
esis underlying the formulation of the DPDE model|18|], and derive the appropriate 



algorithm with arbitrary dependence of the transport coefficients with the temper- 
ature, and energy conservation at every time-step. The Fokker-Planck equation 
describing the evolution of the probability distribution for the ensemble of variables 
describing the state of the system is also derived together with the corresponding 
fluctuation-dissipation theorems. In section III, we obtain the macroscopic equi- 
librium properties of the DPDE model and expressions for the equations of state 
and thermodynamic properties. At the end of this section, we pay attention to a 
particular model, used to perform simulations of equilibrium properties as well as 
heat transport and thermal convection. Finally, in section IV we draw the main 
conclusions from our work. 



2 Dissipative particles with energy conservation 

For our purposes, it is useful to have a physical picture and regard the dissipative 
particles as if they were clusters of true physical particles^ |J, i.e., as particles 
with internal structure bearing some degrees of freedom. Thus, the DPDE model is 
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mesoscopic in nature since it resolves only the overall center-of-mass motion of the 
cluster and ignores the exact internal state of the cluster as a relevant variable. The 
interactions being dissipative and random, however, the total energy of the system is 
not conserved unless the energy exchanged between the resolved degrees of freedom 
and the internal state of the DPD particle is accounted for. We propose here a model 



based on the treatment of thermodynamic and hydrodynamic fluctuations |13|, [T(], 
TTJ, to consistently take into account the energy stored in the internal degrees of 



freedom of each particle, without explicit consideration of any internal Hamiltonian. 



2.1 Definition of the model 

Our model is based on the following assumptions: 

1. The system contains N particles interacting with each other via conservative 
as well as dissipative interactions. The conservative interactions are described 
by the Hamiltonian 

m) = E jf^ + X>(ry) + (2.1) 

where = |r*j — fj\. The Hamiltonian depends on the momenta and the 
positions of all the particles. The particles interact through pair pontentials 
ip(\fi—fj\), depending only on the distance between them, and with an external 
field ip ext (ri) 

2. In addition, the particles can store energy in some internal degrees of freedom. 
The internal energy Ui, with U{ > 0, is introduced as a new relevant coordi- 
nate. The momentum pi, the position r*j together with the internal energy Uj, 
completely specify the state of the dissipative particle at a given instant t. 

3. The particle-particle interaction is pairwise and conserves the total momentum 
and the total energy when the internal energy of the pair is taken into account. 
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4. The internal states of the particle have no dynamics in the sense that they 
are always in equilibrium with themselves. This allows us to define a function 
Sj(ttj). This function can arbitrarily be chosen according to the user's needs, 
except that it is constrained to thermodynamic consistency requirements [D| 



(a) Si must be a differentiable monotonically increasing function of its vari- 
able Ui, so that Ui(si) exists and 6i = dui/dsi exists and is always positive. 

(b) Si is a concave function of its argument. 

Defined in this way, S{ can be viewed as a mesoscopic entropy of the i th particle, 
and 9i can be seen as the particle's temperature. The change in Ui and in Sj 
are related by a Gibbs equation 

9i dsi = diii, which implies 0$ s\ = u\ (2.2) 

where the dot over the variables is used to denote time-differentiation from 
now on. 

5. The irreversible particle-particle interaction is such that the deterministic part 
(in the absence of random forces) must satisfy 

s\ + s) > (2.3) 

where i and j label an arbitrary pair of interacting particles. 

6. Ui, Si and &i must remain unchanged under a Galilean transformation, so that 
these variables are true scalars. 

7. The equilibrium probability distribution for the relevant variables of the sys- 
tem under isothermal conditions is chosen to be 

Pe({Ti}, {Pi}, {Ui}) ~ e -H({n},m)/kT JJ ^{u^/k-^/kT ^ 

i 

where k is Boltzmann's constant and T is the thermodynamic, i.e., macro- 
scopic temperature. The first factor on the right hand side of eq. ( |2.4p con- 
tains the probability distribution for the set of variables {pi}, as given 



by equilibrium statistical mechanics. The second factor on the right hand side 
corresponds, in turn, to the probability distribution for the internal energy 
of the particles as obtained from equilibrium fluctuation theory fT3|j. Effec- 
tively, exp(si(ui) /k — Ui/kT) gives the probability for the i th particle to have 
an internal energy Ui regarding the rest of the system as a heat reservoir. 
The maximum of Si(iii)/k — Ui/kT occurs at 9i = T, in agreement with our 
interpretation of 0j as the particle's temperature. Furthermore, notice that 

where M is the normalisation constant. Once the equilibrium probability 
distribution is obtained, the thermodynamic properties of the model are de- 
termined. 



2.2 Dynamics of the model 

The model defined so far must be completed with a set of equations that explicitly 
establish its dynamical properties. Initially, in view of the pairwise additivity of the 
interactions, we will analyze an arbitrary pair of particles, i and j say, and later on 
give the complete expressions for the iV-particle system. 

The change in position and momentum of the i th particle due to the interaction 
with the j th particle follows from Newton's second law 

h = ^ (2.6) 
m 

ft = pig + FT* + F§ + F§ (2.7) 

where Ff- = —di/;(rij) / dr^ and Ff xt = —dip ext (ri)/dfi are the forces due to the 
conservative interactions. F$ is, by construction, directed along the vector = 
(fj — Ti)/rij and satisfies F$ = —F?. In addition, Ff? stands for the dissipative 
particle-particle interaction force and F^ is the random force associated with the 
former. The total energy of the i th particle is the sum of the kinetic (pf/2m), the 
potential (^{rij) /2 + 'ip ext {f i )) 1 and the internal energy (iii) contributions. Using eq. 
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( 2/7) to compute the change in the kinetic energy, we arrive at the equation for the 



change in the total energy 



r// 2m ' r// V2" ' r V ' J V ' ///' ' V" 'J ' -vj ■ ■>,„■' ' " ' "' 

(2.8) 

The conservation of the total momentum of the pair imposes that its change is 
only due to the action of external force fields, that is, 



? i +? j = Fr t +*F t , (2.9) 

Therefore, we must have that F^' = —F?' R . In addition, conservation of the total 
angular momentum 

fi x £ + fj x = f { x F t ext + fj x Ff- (2.10) 

~* D Ft 

implies that dissipative as well as Brownian forces, F^ ' , must be directed along the 
unit vector fjj. Conservation of the total energy in the particle-particle interaction, 
hi + tj = gives, in turn, the rate of change of the internal energy of the pair 

Ui + *i = ~(pi ~ fi) ■ {FP + F*) (2.11) 

This equation implies that if the total energy is to be conserved, then the change 
in the total internal energy has to be due to the work done by the dissipative and 
the associated random forces. The dissipative forces transfer mechanical energy to 
internal energy, while random forces bring back internal energy to the kinetic and 
potential energies of the particles. From the conservation equations alone, however, 
nothing can be inferred about how the internal energy is distributed among the 
particles so that one has to supply a given model. 

We will assume that the mechanisms driving the change in the internal energy 
of the particles are of two kinds. On the one hand, the work done by dissipative 
and random forces is shared in equal amounts by the particles and is irrespective of 
their temperatures Q{ and 9j. On the other hand, we assume that the particles can 
also vary their internal energy by exchanging internal energy if Qi ^ 9j. The energy 
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tranferred by this mechanism between the particles will be referred to as mesoscopic 
heat flow q®. Associated with this dissipative current, a random heat flow q R must 
also be added. Hence, we write 



Ui = — 

2m 



^-(Pi - Pi) ■ {f% + + «§ + $ ( 2 - 12 ) 



with the requirement q® ,R = —qfi' R to ensure that eq. ( |2.11|) is satisfied. Note that 
the r.h.s. of eq. ( [2.12| ) preserves Galilean- invariance. 

So far, only the conservation equations have been used to find general proper- 
ties to be satisfied by the dissipative forces and mesoscopic heat flows, referred to 
as dissipative currents in a wider sense from now on. In analogy with the Ther- 
modynamics of Irreversible Processes [[Tl]], we can make use of the Gibbs equation, 
eq. (|2.2j ) together with eq. (|2.12|) to find the particle 's entropy production from the 
interaction between pairs, yielding 

* + s > = 1 + % = -Tn ft + Q « "«> ■ ^ + ft " « s (2 ' 13) 

where only the deterministic part of the interactions has explicitly been consid- 
ered, according to point 5. The inequality in this last equation is the expression of 
the Second Law of Thermodynamics, indicating the irreversibility of the particle's 
dissipative interations. From the entropy production equation the so-called thermo- 
dynamic forces can be identified as the factors multiplying the respective dissipative 



currents, and in eq. (|2.13|) . Thus, in the spirit of the Thermodynamics of Ir- 



reversible Processes |T^] we propose a linear relation between the dissipative currents 
and the thermodynamic forces of the form 

p S = 4^ft + £Kv(fi-s) P.14) 

*, . = 4>ft-£) (2.15) 

The functions L$ and L\f are analogous to the so-called Onsager coefficients [fTcR . 
Due to the different tensorial natures of the momentum flux and the heat flux, these 



phenomena are not coupled in eqs. ( |2.14j ) and ( |2.15| ) (Curie's theorem flltl). The 
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mesoscopic Onsager coefficients introduced in eqs. fl2.14|) and (|2.15|) must satisfy 
the following conditions so that the system can reach the proper thermodynamic 
equilibrium 



1. The thermodynamic forces are Galilean- invariant, the dissipative cur- 
rents. This implies that the Onsager coefficients must also be Galilean-invariant. 

2. Microscopic reversibility implies that and L^f must be even functions 
under time-reversal fl4[ . 



3. Since the thermodynamic forces change their sign under the exchange % — > j, 



crucial in the derivation of the transport properties of the system. 



then Lf, and L}y must be invariant under this transformation. This fact is 



In the linear Non-Equilibrium Thermodynamics scheme, L$ and L$ are con- 
stants. However, this choice restricts the expression for the macroscopic transport 
coefficients of the model to a given particular form which, in addition, seldom occurs 



in nature pi] . Thus, we will consider a general dependence of the Onsager coeffi- 
cients in both the distance and the temperatures of the particles 9i and 6j, in 
order to allow the model to describe the temperature dependence of the macroscopic 
transport coefficients. Thus, for convenience, let us rewrite eqs. (|2.14j) and fl2.15|) 
in a different form 



1 

m 

q.. = X^-Oi) (2.17) 



F ij = (ij—hjhj-(pj-Pi) (2.16) 



Written in this way, these expressions are reminiscent of the macroscopic phe- 
nomenological Newton's law for the transport of momentum, and Fourier's law for 



heat transport |10|. Hence, Qj is a mesoscopic analog of the macroscopic viscosity 
and Xij, of the thermal conductivity. Qj and will be referred to as mesoscopic 
dissipative coefficients, from now on. Clearly, the mesoscopic dissipative coefficients 



are related to the mesoscopic Onsager coefficients introduced in eqs. ( [2.14 ) and 
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( |2. 151 ), according to 

L\f = QyCv (2-18) 
4 9) = 9i9j Xij (2.19) 

where O^ 1 = (l/9i + l/9j)/2. In what follows, we will use the mesoscopic dissipative 
coefficients instead of the mesoscopic Onsager coefficients, to base our discussion in 
the most intuitive manner. From eqs. ( 2.18|) and ( 2.19|) one can see that if Qj and 



Xij are constants, then the mesoscopic Onsager coefficients depend on the particles' 
internal energy through the temperatures 9{ and 9j. The presence of temperature- 
dependencies in the mesoscopic Onsager coefficients will introduce subtleties in the 
derivation of the algorithm due to the so-called Ito-Stratonovich dilemma||16|] that 
we will discuss later on. 

The properties of the random terms are also chosen to parallel the theory of 



hydrodynamic fluctuations [111, 1X2 1 . Since and q? are not coupled, we will demand 



that the random terms F$ and be statistically independent. They can be written 
in the form 

'-uVu^jU) and gg = Sign(i-j)A <i Q ii (t) (2.20) 
where the function Sign(i — j) is 1 if i > j and — 1 if i < j, ensuring that q^ = —qfy 



The scalar random variables and are stationary, Gaussian and white [TJ, 15 
with zero mean and correlations 

(Fij(t)Fki(t')) = (QiMQki(t')) = (6 lk 5 jt + 8 a 5 jk ) 6(t - t') (2.21) 

(^i(0S«(O) = (Qij{t)Fki{t')) =0 (2.22) 

Tij and Ay are functions to be determined later. Note that Qj and are, respec- 
tively, 7c<Jd and au R in ref.|J. 



2.2.1 Derivation of the algorithm 



As we have seen, the dynamics of the model is described by the set of Langevin 
equations, eqs. Q2.6Q , Q2.71) and (|2.12|) , the definition of the dissipative currents, 
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eqs. (|2. 16 ) and Q2.17 ), and the random forces, eqs. ( |2.20| ). To derive an algorithm 
describing such a stochastic process, let us integrate these Langevin equations for a 
small increment of time St and retain terms up to 0(5t). One obtains 



r' 



Pi 



U; 



n H bt 



rn 



(2.23) 
(2.24) 



■it; 



E 



[(Pi - Pi) ■ rij} 2 + \j (% 



OA \ St 



+ E { 2^-03 - ) • V* ^ 1/2 ^S? + sign(* - i)A^. <^/ 2 n 



(2.25) 



where r/ = fi(t + 5t), p/ = pi(t + St), and = + 5t) while r i; ^ and are 
the value of these functions at time t, and r'„- and A'.- stand for the value of these 
functions when their arguments are calculated at the time t + St. The integrals over 



dr r y (T)r y [T].Ftf (t) = h 3 {t)T l3 [t + 5t] c^ 1 / 2 ng> 



the random terms have been interpreted as [16] 

rt+St 

Jt 

in eq. ( ggg ), and 



;2.26) 



JT" rfr (^(r) - fl(r)) • f^r) r y [r] ^(r) 

(p,-(t + 5t) - Pi(t + <Jt)) • r„[t + ft] ft 1/2 ftjf (2.27) 
rfr Sign(* - j)Ay[r] Qy(r) = Sign(* - j)K l3 [t + ft] ft 1 / 2 (2.28) 



respectively, in eq. ( 2.25|) . In these equations, we have defined the random numbers 

1 rt+St , . I rt+St 

v* Jt dr Tii ^ and s Jl dr 



lJ ft 



(2.29) 



which are Gaussian, with zero mean and correlations (f^fij^) = (fiy fij^ ) = 
(Mii + <y«<5>-fc) and (fi^fifc?) = 0, according to eqs. (gjjD and (g^2|) . Eqs. ( gjg ), 
( p.27|) and ( |2.28| ) express the fact that the amplitudes of the random forces may 
depend on the fast variables themselves (ui through the temperature 0$, in this 
case), introducing the so-called Ito-Stratonovich dilemma ||16|]. Here, we have chosen 
an interpretation rule that is neither that of Ito nor Stratonovich, but has the 
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property that the resulting stochastic process, once the amplitudes of the random 
terms have been fixed, satisfies detailed balance up to 0(5£)p7|. Note that in eqs. 
Q2.26 ), ( 2.27 ) and ( |2.28 ) the slow variable is taken at t instead of at t + St, since 
the difference between both considerations is of order 0(5t 3 ^ 2 ); only fast variables, 
whose increments contain terms of the order 0(St l l 2 ), need to be taken at the latest 
time since they introduce terms of order 0(St) in the algorithm. That the stochastic 
process introduced by eqs.( |2.23| ), ( |2.24| ), and ( |2.25| ) satisfies detailed balance can be 
verified from the functional form of the resulting Fokker-Planck equation, since this is 
a property of the transition probabilities and the equilibrium distribution|16[. Other 
algorithms have been proposed that explicitly incorporate time-reversibility in the 
equations of motion for the particles ||24|| , which eventually leads to the resulting 



stochastic process satisfying detailed balance^]. The forward interpretation rule, 
however, permits us to straightforwardly relate the differential equations for the 
change in the particle's variables (Langevin equations) with an Euler algorithm that 
will lead to the proper thermodynamic equilibrium for the system. 



To obtain a causal or explicit form of the algorithm defined in eqs. ( |2.23| ), ( |2.24| ) 
and (|2.25|) , we expand the right hand side of these equations in powers of St and 
retain terms of up to first order. We then obtain 



5fi = — St 
m 



Spi 




d d 



'3 , 



13 ' \ m 2m 13 \ dui ' dun ' 13 13 



+ f l3 Sign(z - j) A,, (£-- S-\ TM^nW 



13 \d Ul du 3 ) 13 13 13 



1 fCi 



1 



d d 



~f. 1 2m \ m 2m \ dui duj 

1 

2m 



! —(^-^.f^^fSign^-j; 



(2.30) 



St + ^f^St^n^ (2 .3i) 

3+i 



, d 9 \ » a f 9 d , 

r« l — + — I A** - A iA I - — — I r 



1 dui ' duj J 1 ^ 3 %3 \ duj dui ' * ' ' 



A; 



d d 



2 o(P) 2 



m 



13 \ duj dui J %3 13 
+ E { 2^ (PJ - Pi) ■ raVii st 1 ' 2 njj> + A l3 st" 2 ng> } 



(2.32) 
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Eqs. ( |2.30|) , (|2.31| ) and (|2.32|) is an explicit algorithm, equivalent to the implicit 
algorithm given in eqs. Q2.23j ), ( |2.24| ) and ( |2.25| ) up to order 0(5t). By using 
standard methods[16|, from eqs. Q2.3C ), ( 2.31| ) and (|2.32 ) one obtains the Fokker- 
Planck equation 

^ p ({ri\, m, M) = L c ° n p({n}, m, { Ul }) + L^Piin}, m, k» (2.33) 

where the convective operator, L con , and the diffusive operator, L dl f , are defined as 





+ 



N 

-£ 

1=1 

d 



Out 

N 

L di f = J- 



A . B. + — . p ext 

dfi m dpi 1 

Cij 



N 



F% + 7T'"'j'"'J ■ (Pj ~ Pi) 



m 



+ 



2m 2 
d 1 



'j-Pi)-hj\ +\ ij (9 j -6 i ) 



_d_ 

dui 



l^Pi - Pi) ■ ^ •/■ J '->j ■ Ai + 2 ij [ ^7 



where we have in addition defined the operator 

Note that the algorithm giving the aforementioned Fokker-Planck equation is not 
unique. Effectively, the terms of order O(St) in eqs. ( |2.30| ), ( |2.31| ) and ( |2.32| ) can be 
replaced by their averages over the random terms, giving a much simpler algorithm 



satisfying the same Fokker-Planck equation, eq. (|2.33| ). Such a procedure was 



adopted in the first version of the DPDE model [18]. However, replacing the drift 
terms by their averages means that the energy is only conserved in the mean but not 
at every time-step. The implicit algorithm shown in eqs. ( |2.23|) , ( |2.24j) and ( |2.25| ) 



is thus a more compact form satisfying the desired properties of the DPDE model 
described. 

The amplitudes of the random terms are derived by imposing that the equilibrium 
distribution function given in eq. (|2.4j ) is a stationary solution of eq. (|2.33|) , thus 
yielding the fluctuation-dissipation theorems 

r| = 2fc4 p) =2fcCi 3 -8 ij (2.37) 
Aj = 2k L[f = 2k$i8j Xij (2.38) 
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Eq. (|2.33|) , in view of the definitions ( |2.34|) , (|2.35|) and (|2.36| ), together with the 



expressions for the amplitudes of the random terms ( |2.37| ) and ( p.38| ), satisfies the 



expected detailed balance condition [14 



The fluctuation-dissipation theorems given in eqs. (|2.37p and (|2.38|) are the gen- 



eralization of the results previously obtained in ref . |T8| , to the case of an arbitrary 
temperature-dependence of the mesoscopic coefficients. We have to stress the fact 
that the amplitude of the random force depends only on particle variables and 
that it is independent of the macroscopic state of the system (the thermodynamic 
temperature, for instance). This point is a crucial difference between the DPDE 
model and the previous isothermal DPD models, with no energy conservation. For 
the isothermal DPD model, the amplitude of the random force is proportional to 
the thermodynamic temperature of the system T which must be previously speci- 
fied. This particular feature of the actual DPDE allows one to model temperature 
gradients in the systems as well as heat flows. Furthermore, since the dynamics is 
based on particle properties, DPDE can deal with systems either in contact with a 
heat reservoir or with isolated systems. 

Having set the fluctuation-dissipation theorems that fix the amplitude of the 
random terms, the algorithm describing the DPDE process is completely specified. 



3 Thermodynamics of the DPDE model 

The properties of the DPDE system defined so far are such that the system evolves 
irreversibly towards a final equilibrium state. If the system is in contact with a heat 
reservoir, such an equilibrium state is characterized by the probability distribution 
given in eq. ( [2.4|) . Since this fact allows us to define a thermodynamics for the system 
of dissipative particles, we will devote this section to the calculation of general results 
valid for any model. 
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Let us introduce a partition function in a classical sense, Q(T, V, N), as 

Q(T,V,N) = j {Xldvidndu^ e-^WM/M-JJe^Vifc-m/w 

= Q H (T,V,N)Q int (T,N) (3.1) 

where Q# refers to the partition function of the " hamiltonian" part of the interaction 
and Qint, to the partition function of the internal degrees of freedom. The factor 
l/h 3N has been introduced in analogy with the partition function of a physical 
system. The equivalent Free Energy for the system of dissipative particles can thus 
be obtained as a sum of two contributions 

F{T, V, N) = -kT In Q H (T, V, N) - kT In Q mt (T, N) (3.2) 

The former is related to the hamiltonian part of the dynamics and the latter, related 
to the internal energy dynamics. Qh, can be obtained from the Hamiltonian of the 
system given in eq. (|2.1|). One gets 

Q„(T,V,N) = J (jldRdr^ e -»<«>.«»/«' 

1 f nmt -^^^ . ^Z H (T,V,N) 



A 3N N \ J 11 1 A^NV 

(3.3) 

where A = hj 1 \phxmkT is the de Broglie wavelength and Zh(T, V, N) is the so-called 
configurational integral^T7\ . The internal partition function can also be worked out. 
Effectively, one has 

N 



Q int (T, N) = jY[duill e^/k-vi/kT = J 



due s(u)/k-u/kT 



[Qi(T)f (3.4) 



hence factorizing in one-particle partition functions, depending only on the temper- 
ature. 

Using these expressions, the free energy takes the form 

F(T, V, N) = kTN{ln p + 3 In A - 1) - kT In -^^^ - kTN In Q x (T) (3.5) 

Note that Qm is a function of the temperature and the number of particles, but 
is independent of the volume. Therefore, the pressure of the system is completely 
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determined by the Hamiltonian interactions between particles, as one could have 
intuitively guessed. The internal energy, however, contains contributions from all 
the terms. 



The explicit knowledge of the equations of state of the system of DPD particles 
are very useful, since they can be used to compare with the equations of state of 
a real system when simulating its behaviour and properties in a computer. In the 
rest of this section we will derive general equations of state for the pressure and 
also for the internal energy of the system, in terms of microscopic parameters of 
the model. The analysis given here will follow that of ref. flTj, where the details of 
the calculation can be found. Let us introduce the one-particle and the two-particle 
distribution functions, according to the relationships 

p(rO = ^E^-n))> ( 3 - 6 ) 

for the former, and 

p( 2 )(f,f/) = (£5(r- fi)5{f/ - fj)\ (3.7) 

for the latter. The averages refer to the equilibrium probability distribution for 
the DPD system given in eq. fl2.4|) . For a homogeneous and isotropic system, 
these expressions can be further simplified. Effectively, the one-particle distribution 
function is the mean density p(r) = N/V and the two-particle distribution function 
reduces to 

pW(\r-rf\) = p 2 g(\f-rf\) (3.8) 
where this last equation is in fact a definition of the pair-distribution function g(r). 

The pressure equation can be obtained from the so-called virial equation 

^i + 4(£-4s/ M ) (3 - 9) 



where j3 = 1/kT. After some algebra, one arrives at the final expression fTF 



p 3 Jo dr 
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From this expression one can see that the actual model of DPDE particles can 
exhibit the same phase behaviour as a fluid of pair interaction potentials, thus gas, 
liquid and solid phases can be modelled. The actual DPDE shares this property 
with the older versions of the DPD model. 

The internal energy of the system can be obtained by performing the equilibrium 
average of all the contributions to this quantity 

U = J df(j2 - + / drdr/ ^ HnjW- ?i)5(r/ - 

+ J drl^u^r-r^ (3.11) 

Again, after some algebra, we obtain 

U = ^NkT + + Ne u (3. 12) 

where the first term on the right hand side of this equation is the ideal gas contribu- 
tion. The second term is the interaction potential energy contribution, and is given 
by 

N r 

# = 2vr— J drr 2 g(r)i)(r) (3.13) 

The third term is the contribution due to the internal energy of the particles to the 
macroscopic internal energy of the system. It is evaluated from the integral 

f^°duue s ^/ k - u / kT 

Note that the thermal capacity of the system is obtained by differentiation, dU/dT) v . 
This can be used to relate the function s(u), that has to be supplied to the model, 
with the thermal behaviour of a given real system. These expressions cannot be fur- 
ther developed without explicit calculation of the pair distribution function. How- 
ever, they can shed some light on the relevant aspects of the thermodynamic equi- 
librium in the DPDE system. 

To end this section, let us introduce a non-equilibrium extension of the free 
energy, given by the expression 

F[P] = J dX P(X,t){f(X) + kT\nP(X,t)} (3.15) 
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which is defined as a functional of the actual probability distribution P, and where 
X is a point in phase space, {{Pi}, {^}, {ui}), and the function f(X) is defined as 

f(X) = Hdpi}, {fi}) T Si (ui)) = -kT In P e (X) (3.16) 

i 

Using this expression for the function f{X), the eq. fl3.15|) can be cast in a more 
convenient form 

P(X,t) 



T[P] = J dXkTP{X,t) In 



Pe{X) 



(3.17) 



This functional T satisfies an if -theorem, dT jdt < 0, as follows from differentiation 
of eq. ( |3.17| ) with respect to time, with the use of the Fokker-Planck equation given 
in eq. (|2.33|) together with the fluctuation-dissipation theorems (|2.37|) and (|2.38|). 



3.1 Analysis of a particular system 

In the remainder of this section we will analyze a particular model and show some 
simulation results. 

Let us consider the case in which the DPDE particles have no interaction poten- 
tial forces. The definition of the model requires a particle equation of state relating 
the particle's entropy with the particle's internal energy, s(ui). The simplest model 
is to assume that the particle temperature is a linear function of the particle internal 
energy of the form Ui = 4>6i, where is a constant playing the role of a particle's heat 
capacity. We will take to be larger than the Boltzmann constant fc |p3[| . Hence, 
the particle internal energy probability distribution takes the form 

P e ( Ui ) ~ e^/k-^/kT = u 4,/k e - Ui /kT ( 3 _ 18 ) 

In the absence of an interparticle interaction potential the position probability dis- 
tribution is uniform, and the momentum distribution is that of Maxwell-Boltzmann 

p 2 

P e (Pi)~e-^ (3.19) 

The model is then analogous to a general ideal gas and, hence, it corresponds to a 
single gas phase, as seen from the resulting pressure equation 

p = kTp, (3.20) 
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in view of eq. (|3.10|) . The second equation of state determining the thermodynamic 



state of the system, follows from eq. G3.12| ) and gives the variation of the macroscopic 



internal energy with respect to the state variables. For our particular model, one 

gets 

U = ^NkT + NkT ( | + lj (3.21) 

In this equation, the first term is the ideal gas contribution to the macroscopic 
internal energy, due to the translational degrees of freedom of the particles. The 
second term, 0T, is the internal energy stored per particle. Note, however, the 
additonal kT on the right hand side of eq. ( |3.21| ) , that comes from the extra degree 



of freedom due to the fluctuations in the internal energy of the particles. The 
heat capacity can be calculated by differentiating eq. (3-21[) with respect to the 
temperature, giving 

C v = NkU + ^ (3.22) 



Let us point out that eq. ( 3.22 ) is a relationship between macroscopically measurable 



magnitudes of any real system, C v , with model parameters </>. 

As far as the dynamic properties are concerned, the particle friction and thermal 
conductivity are given by the expressions 

Qj = Co (l - for < r c (3.23) 

where Co an d are constants giving the magnitude of the mesoscopic friction and 
thermal conductivity, respectively, and and r\ are the respective ranges of the 
dissipative interactions. The functions Qj and Ay vanish if > and > r\, 
respectively. Although other choices could be made for the spatial dependence of 
Qj and we have made here the usual choice Jl|. 

The explicit derivation of the macroscopic transport properties of the system will 
be treated elsewhere. Here, however, we give results that can be easily calculated 
by considering that the transport of momenta and energy is dominated by the dis- 
sipative interactions. Such a limit must be reached under conditions of high density 
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system or by choosing large particle friction and thermal conductivity parameters. 
By considering the particles as frozen and analyzing the transport of mesoscopic heat 



through a hypothetical plane dividing the system into two parts[^6|, one obtains for 
the macroscopic thermal conductivity in these limiting conditions 

n 2% m , . / r ..\ 2 

T\. 



A 



^r^^K 1 -?) 2 ^ (3 - 25) 



A similar calculation for the transport of momentum yields the shear viscosity 



2 



^^/Vr 4 C (l-^j g{r) (3.26) 

which coincides with the result obtained in ref.|| when the particle friction coeffi- 
cient is considered as very large. Performing the integrals given in eqs. ( |3.25 ) and 



( |3.26| ) by using the fact that in equilibrium and in the absence of pair interaction 
potentials, g(r) = 1, we then obtain for the viscosity 



= ^f-rko (3.27) 
' 1575 CS 



For the macroscopic thermal conductivity, one gets 

A = —^rl (3.28) 
315 T 2 v ; 

The analysis of the simulation results with steady state heat conduction shows agree- 
ment with the functional dependence given in eq. fl3.28|) . A future publication will 



be devoted to a deeper analysis of the transport properties of the DPDE model. 

The velocity of sound is theoretically obtained by means of a thermodynamic 
calculation 

dp) s C v dp) T (3.29) 
For the case under discussion, one obtains for the speed of sound the ideal gas result 

+ 7fe/2 kT 
° ~ + 5fc/2 m [ } 



In order to carry out the simulations, we have introduced dimensionless variables 
suitable for a proper interpretation of the results. We have defined a temperature 
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of reference Tr to be used as the scale of temperature. Thus, T = TrT* and 
Qi = Tr6**, where the asterisk is used to denote a dimensionless variable from now 
on. The momentum of the particles is made dimensionless according with pi = 
p^y/2mkTji, using the characteristic value of the momentum in a system with no 



externally imposed flow and thermal fluctuations: J (pf). The penetration depth for 



the momentum defines a characteristic length scale / = y2mA?T#/£o- Thus, r*j = r*l. 
The characteristic scale of time is the relaxation time for the particle's momentum, 
that is t = t*m/(o. The characteristic scale for the particle's internal energy is 
chosen to be 0Tr so that Ui = w*0Tr. When the algorithm is written in terms of 
these dimensionless variables it can be seen that there are only two independent 
dimensionless parameters in the model described in this section, that is 

B = \ (3.31) 
C = ^ (3.32) 

Com 

B measures the relative magnitude of the fluctuations in the particle's internal en- 
ergy with the characteristic particle energy. It is thus convenient that B be small. C 
is the ratio between a relaxation time for momentum decay m/( , and the relaxation 
time for the decay of a particle's internal energy fluctuations (pT^/L q \ Note that, 
to avoid negative values of the particle's internal energy during the integration of 
the algorithm (eqs. ( gjg ), (gjg and (Hj)), in the presence of externally imposed 
temperature gradients, one can roughly estimate that 

C ( Y r f) P* (r*\V*T*) St* < 1 (3.33) 

where p* is the dimensionless particle number density and St* is the dimensionless 
integration time step. The fluctuating part of the energy can also lead to negative 
energy values if the condition 

47T 



VBCSt* ( — rf ) p* < 1 (3.34) 

is not satisfied. 



Two aspects of the DPDE model have been investigated using simulations in 
three dimensions. In the first place, the thermodynamic consistency has been 
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checked by measuring the equilibrium distributions for a system in contact with 
a heat reservoir. Figs. 1 and 2 show the momentum and particle's internal energy 
distributions of a system of N = 30000 particles in contact with two walls at a 
temperature T* = 1.1, while periodic boundary conditions are chosen for the other 
two dimensions in space. The factor B has been set equal to 0.01, C — 17 and the 
cutoff lengths and are both set equal to 1.24. The size of the system is chosen 
such that p* = 1, so that the dimensionless lateral size of the box is L* = N 1 ^ 3 . 
The number of particles interacting simultaneously with a given particle is thus de- 
termined by the size of r* x and rJ. For our values, we can roughly estimate that 8 
particles interact at one time with a given particle. 

The simulation has been initialised by a random distribution of particles at rest 
and a dimensionless particle temperature 9* = 1. At the initial stages of the simula- 
tion, the system is cooled down to a temperature of about 0.98 by the transformation 
of internal energy into kinetic energy due to the action of the random forces, in a 
time scale t* of about 1. The difference in temperature between the wall and the bulk 
provokes a heat flow from the walls to the particles (see Fig 3). This process takes 
place in a much longer time scale, found of the order of 100, although this is clearly 
related to the overall size of the system. In Fig. 3 we show the time evolution of 
the mean temperature of the particles. After equilibration, this mean temperature 
reaches the same temperature as that of the walls. In addition, the three transla- 
tional degrees of freedom are found to satisfy the probability distribution given in 
eq. ( |3.18| ) which, expressed in dimensionless form reads 

Peivla) = 7^ e " P! * i/T * ( 3 - 35 ) 

where p* a stands for any of the three components, a, of the momentum of the i th 
particle, and T* is equal to the temperature of the wall, 1.1. In Fig.l, simulation 
and theoretical results are shown together. We find excellent agreement between 
simulated and theoretical distributions which, in addition, are rather insensitive to 
the time step, provided it is small (5t* = 0.01 in this simulation). In equilibrium, 
the average internal energy of the system is found to be U* = 33828 ± 1. The same 
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property calculated by eq. (|3.21| ) (written in dimensionless form) yields the value 



U* = 33825, showing excellent agreement with the simulation value. 

The degree of freedom represented by the particle's internal energy also satisfies 
the probability distribution given in eq. (|3.19|) . This distribution can also be cast 
in dimensionless form, yielding 

*1/B _,,*/BT* 

where here T stands for the Euler Gamma function, and T* is the temperature of 
the walls. Fig. 2, shows the theoretical predictions and simulation results for this 
probability distribution. Once more, the agreement is excellent. This result is non 
trivial in view of the fact that the functional form of the probability distribution for 
the energy is arbitrary, and depends on the choice made for the particle equation 
of state, Ui = (j)9i in our case. We have found, however, that this probability 
distribution is much more sensitive to the time-step (5t* = 0.0001 in the simulation 
shown in Fig. 2). For a time step 5t* = 0.01, the deviations are less than 1%. 

Finally, the pressure in the simulation can be obtained from the average force 
that the particles exert on the walls. For the previously mentioned system we have 
obtained the pressure from the simulation and found a value of P* = 1.137 ± 0.005. 
The theoretical value given from the equation of state eq. ( |3.20|) is P* = 1.1. Again 
we find a good agreement between simulation and theory. 

Simulations of closed systems (with periodic boundary conditions in the three di- 
mensions) have also been performed (this is possible because the algorithm depends 
on particle properties only). The DPDE system tends in this case towards a final 
thermodynamic equilibrium satisfying the aforementioned probability distributions 
for the one-particle variables. In this case, however, the probability distribution 
of the complete system is not separable as it was in the previous case. The final 
temperature attained for the system is in excellent agreement with that predicted 
from eq. (|3.21|) , provided that the total energy of the system is known. 



The second aspect analysed is the ability of the model to represent non-equilibrium 
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features of fluids under temperature gradients, an unattainable problem for the 
older DPD algorithms. We have performed simulations of a DPDE system with 
N = 100000 particles in a cubic box of lateral size iV 1//3 , with four walls and pe- 
riodic boundary conditions in the remaining direction. The walls were considered 
to be a two-dimensional surface which exert a repulsive force in the orthogonal di- 
rection on the particles. The range of this force is of the order of the range of the 
particle-particle interaction. As far as the wall-particle dissipative interactions are 
concerned, we have considered as being analogous to the particle-particle interac- 
tions. Hence when a particle is within the interaction range of a wall, it exchanges 
heat and exerts a force orthogonal to the wall as if the latter were a particle of in- 
finite mass and heat capacity. Therefore, the walls act as heat reservoirs at a given 
specified temperature, and non-stick boundary conditions apply. 

In the simulations, the two walls in the x-axis were kept at fixed temperatures of 
T£ = 2.8 and T* = 0.8, while the other two walls were adiabatic and were placed in 
the z-axis, along which a gravity field was imposed on the particles. Qualitatively, 
the gravity was weak enough not to cause an excesive density gradient along the z- 
axis, and the particle's thermal condutivity was chosen to be very small to emphasize 
the convective effects. Initially at rest, the system spontaneously evolved towards 
a stationary convection roll as seen in Fig. 4, where the arrows stand for the fluid 
velocity field. The velocity field has been calculated by dividing the space into 
boxes and averaging the velocity of the particles inside the box at a given instance 
of time. The fields obtained in this way show a rather marked variability due to 
the inherent stochastic nature of the algorithm and the limited number of particles 
inside a given box at a given time. The results shown in Fig. 4 are thus the 
result of averaging in time as well as in the ^/-direction, to have better statistics. 
In addition, the temperature profile has been obtained. We see in Fig. 5 that a 
non-linear temperature profile, due to the compressibility of the model, exists and is 
slightly tilted due to the convective flow. These results are in qualitative agreement 
with convective rolls observed in boxes and in numerical solutions of the Navier- 
Stokes equations for incompressible systems under the conditions described here. 
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The estimated Rayleigh number is of the order of 1000, by comparison with the 
numerical solution of the Navier-Stokes equations. In Fig. 6 we show the density 
profile given by the model, in which one can appreciate the combined effects of the 
temperature gradient and gravity field, merging into a diagonal density gradient, 
due to the compressibility of the system. It is interesting to note how this density 
profile affects the convection pattern in the system (Fig. 4), causing a displacement 
of the vortex centre towards the denser region, and also inducing a faster motion in 
the upper layers, as compared with the bottom of the box. 



4 Conclusions 

Various aspects of the Dissipative Particle Dynamics with Energy Conservation algo- 



rithm, already introduced in ref.JTBJ, have been treated here in depth. In particular, 
emphasis has been placed on two major points. Firstly, the original DPDE algorithm 
has been extended in two ways: to incorporate arbitrary temperature-dependencies 
in the transport coefficients, and to assure that the algorithm conserves the energy 
at every time-step rather than in the mean, as in our previous algorithm[ll8|. Sec- 
ondly, in this paper we have studied the thermodynamic properties modelled by the 
DPDE algorithm such as the free energy and equations of state, and shown agree- 
ment between the simulated and theoretically predicted probability distributions. 
To demonstrate the ability of the DPDE algorithm in simulating the hydrodynamic 
behaviour of fluids under non-equilibrium conditions, simulations have been carried 
out for a system with a temperature gradient orthogonal to a gravity field. The 
results show the expected convective pattern for the velocity field, as well as the 
corresponding tilted temperature profile. These points stress the inherent features 
of the new DPDE algorithm as compared with the isothermal versions. 

With respect to the first major point, and with a view to the use of the DPD 
methodology for the simulation of the dynamics of real fluids, it was necessary to 
incorporate arbitrary dependencies of the transport coefficients in the temperature. 
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We have constructed an algorithm where the mesoscopic friction and heat conduction 
coefficients, Qj and Ajj, respectively, can depend on the particle's temperature. One 



can compare eq. (|3.23|) with eq. ( |3.26|) , to see that both coefficients are independent 
of the temperature, neither the particle's nor the thermodynamic temperatures. 
Contrarily, the functional form of eq. (|3.24|) gives rise to the 1/T 2 dependence of 
the thermal conductivity coefficient. Of course, since eqs. ( |3.25|) and (|3.26|) are 
approximate, there must exist additional temperature dependencies that may be 
due either to g(r) or to the kinetic transport of momentum and energy. Note that, 
the latter has been neglected in our derivation of the transport coefficients. 

Another aspect worth mentioning in this context is our derivation of an implicit 
algorithm. This derivation has two important properties. On the one hand, the al- 
gorithm can be directly obtained from the Langevin equations used to formulate the 
problem of the DPD particle dynamics. We have simply introduced an interpretation 
rule of the random terms, when integrating the equations of motion in a 5t, which 
is different from the usual Ito-Stratonovich interpretations. On the other hand, 
the resulting algorithm, with the proper fluctuation-dissipation theorems, straight- 
forwardly satisfies detailed balance which is required for the model to behave in a 
thermodynamically consistent manner. Furthermore, since this algorithm has been 
directly obtained from the Langevin equations (|2.6|), ( |2.7| ), and Q2.12| ), it satisfies 



energy conservation at every time step, as do the Langevin equations themselves. 

We have also indicated that a family of different algorithms can lead to the same 
Fokker-Planck equation. Since the process is Gaussian, it is completely determined 
by the first and second moments of the probability distributions for the random 
variables. Thus, the defining trends of this family of algorithms is that these first 
and second moments be the same for all of them, up to the order of validity of the 
algorithm itself, 0(St). As an example, let us average eq. ( |2.31| ) with respect to the 
random number fij? , giving 



^ + (| + i r «fe + ^-l r «)(ft-p.)-^ 



The same result is obtained from our previous algorithm proposed in eq. (16) of 



(4.1) 
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ref. ||18|| . This demonstrates that indeed the firsts moments of the variation of the 
momentum for both algorithms are identical. Of course, the same procedure can be 
repeated for all the first and second moments of all the increments of the variables 
in a time-step. Thus, despite the fact that the macroscopic properties of both algo- 
rithms are the same, it should be emphasised that the present algorithm preserves 
energy conservation at every time-step while all the other possible algorithms only 
do so in the mean. This can be verified by summing all the contributions to the 
energy for a pair of particles in a time-step. This point is of crucial importance 
when analysing the dynamic behaviour of the system, in particular when the exact 
conservation of the energy is checked during the simulations. 

Regarding the second major point, we have emphasized in this article the ex- 
isting tight relation between the overall properties of the DPDE model developed 
here and the macroscopic properties of real systems. This tight relationship will 
allow DPDE models to simulate the dynamic behaviour of real fluids, with respect 
to the equilibrium and transport properties of both momentum and energy. When 
the energy conservation is introduced, the existence of a Thermodynamic behaviour 
of the DPDE system becomes apparent. We have seen that the hypothesis leading 
to the present algorithm can reproduce not only the Maxwell-Boltzmann probability 
distribution for the momenta, but also the arbitrary probability distribution for the 
particle's internal energy. Furthermore, we have introduced a free energy from a 
partition function related to the system probability distribution. In the calculation 
of this partition function, the particle's internal energy Ui has been considered as 
an additional microscopic degree of freedom. From this formulation, we can obtain 
equations of state in terms of the parameters of the model. In particular, we have 
obtained a relationship between the internal energy of the whole system and its 
macroscopic temperature. This equation of state is useful since it allows one to re- 
late the total energy introduced in a closed DPDE system with the final equilibrium 
temperature. The agreement between the predicted temperature and the simulated 
one is excellent. The pressure equation for this system has also been obtained, 
showing once more excellent agreement between the simulated and theoretical pre- 
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dictions. The pressure in a DPD system in general is constrained by the particle 
number density, which is much smaller than that of the corresponding real system 
being modelled. This is due to the implicit coarse-graining in the model, in which a 
packet of physical molecules is represented by a single DPD particle. Clearly, this is 
an area that needs to be further studied. However, for low Mach number flows, the 
DPDE fluid can be regarded as effectively incompressible, and thus a Boussinesq-like 
approximation is of use. Under these conditions, the thermal dilatation coefficient 
is the one that has to be properly modelled. In the particular system analyzed here, 
for instance, this coefficient, a, is that of an ideal gas 



With a proper choice of the interaction potential this coefficient could be tuned ac- 
cording to the kind of fluid to be modelled. Note the excellent qualitative agreement 
found between the non-equilibrium simulations of a convective flow presented here 
and real fluid flows under equivalent conditions. These simulations can be theo- 
retically analysed in light of the Boussinesq approximation. Finally, we have also 
identified a non-equilibrium free energy which acts as a Lyapunov functional for the 
dynamics of the system. Therefore, as seen in the simulations the DPDE system has 
an inherent irreversible tendency towards a thermal equilibrium when no external 
forcing exists. Moreover, other quantities have been theoretically derived such as 
the viscosity, the thermal conductivity, and the speed of sound. An in depth analysis 
of the derivation of these quantities and the consequences that can be drawn from 
its dynamic behaviour lie beyond the scope of the present work. 

Therefore, the inclusion of the internal energy in the DPD algorithm reinforces 
the internal consistency of the model and opens interesting perspectives aimed at 
the use of these models in the simulation of the macroscopic-mesoscopic behaviour 
of simple and complex fluids. 




(4.2) 
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Figure Captions 



Fig. 1: Equilibrium momentum probability distribution expressed in reduced vari- 
ables, under the conditions (T* = 1.1, N = 30000, C = 17, B = 0M,r x = r c = 
1.24). The solid line stands for the theoretical prediction given in eq. (|3.35| ) and 
the dots are simulation results. The time-step is St* = 0.01. 

Fig. 2: Equilibrium particle's internal energy distribution expressed in reduced 
variables, under the conditions (T* = 1.1, N = 30000, C = 17, B = 0.01,r A = r c = 
1.24). The solid line is a plot of the theoretical distribution given in eq. ( 3.36Q . The 
time-step used in this simulation is 5t* = 0.0001. 

Fig. 3: Equilibration of the temperature starting from a random configuration under 
the same conditions as given in Fig. 1. 

Fig. 4: Velocity field for a system under a temperature gradient orthogonal to 
a gravity field. The wall temperatures were chosen to be = 2.8, on the left 
hand side, and T* = 0.8 on the right. The gravity field was g* = 0.01, directed 
downwards. The system contains N = 100000 particles in a cubic box of side 
iV 1 / 3 , with adiabatic walls in the top and the bottom sides of the cube and periodic 
boundary conditions in the remaining direction. The additional parameters were 
chosen B = 0.1, C = 0.001, r c = r A = 2, with a time-step 6t* = 0.01. 

Fig. 5: Isothermal lines for the system described in Fig. 4. The numbers stand for 
the thermodynamic temperature T* for a given line. 

Fig. 6: Constant density lines for the system described in Fig. 3. The numbers 
represent the values of the density p* for a given line. The distortion near the walls 
are caused by repulsive potentials exerted by the walls used to confine the system. 
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